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Algorithms for the computation of the real zeros of hypergeometric func- 
tions which are solutions of second order ODEs are described. The algo- 
rithms are based on global fixed point iterations which apply to families of 
functions satisfying first order linear difference differential equations with 
continuous coefficients. In order to compute the zeros of arbitrary solutions 
of the hypergeometric equations, we have at our disposal several different 
sets of difference differential equations (DDE) . We analyze the behavior of 
these different sets regarding the rate of convergence of the associated fixed 
point iteration. It is shown how combinations of different sets of DDEs, 
depending on the range of parameters and the dependent variable, is able 
to produce efficient methods for the computation of zeros with a fairly 
uniform convergence rate for each zero. 

AMS subject classification: 33Cxx, 65H05 
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1. INTRODUCTION 

The zeros of hypergeometric functions are quantities which appear in a vast num- 
ber of physical and mathematical applications. For example, the zeros of classical 
orthogonal polynomials (OP) are the nodes of Gaussian quadrature; classical OP 
(Hermite, Laguerre and Jacobi polynomials) are particular cases of hypergeomet- 
ric functions. Also, the zeros of Bessel functions and their derivatives appear in 
many physical applications and there exists a variety of methods of software for 
computing these zeros. 

However, an efficient algorithm which can be applied to the computation of all 
the zeros of any hypergeometric function in any real interval (not containing a 
singular point of the defining ODE) is still missing. 
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In [7, 1] methods were introduced which are capable of performing this task for 
hypergeometric functions which are solutions of a second order ODE; an explicit 
Maple algorithm was presented in [4]. The starting point of the methods is the 
construction of a first order system of differential equations 

y'{x) = a{x)y{x) + 5{x)w{x) 
■w'{x) = I3(x)w(x) + ^{x)y(x), 

with continuous coefficients a(x), P{x), 7(x) and ^(a;) in the interval of interest, 
relating our problem function y(,;:) with a contrast function whose zeros arc 

interlaced with those of y{x). Typically, the contrast function w{x) satisfies a second 
order ODE similar to the second order ODE satisfied by the problem function. 

Given a hypergeometric function y{x) there are several known options to choose 
as contrast functions w{x). As an example, considering a Jacobi polynomial 

y{x) = PjT'^Hx) = + 2Fi(-n, n + a + P + 1; a + 1; {1 - x)/2) (2) 

we could take as contrast function Wov{x) = P^'^i\x) but also Wd{x) = -^Pn"'^\x) 
is a possible choice. Both contrast functions are 2Fi(a, b; c; x) hypergeometric func- 
tions with parameters a, b, c differing by integer numbers from the parameters of 
the problem function y{x). 

When the contrast function in the previous example is Wop{x) = P^"'f^(a;) the 
first order differential system is related to the three term recurrence relation for 
Jacobi polynomials. It may seem that this is a natural differential system to con- 
sider. However, it was numerically observed that the fixed point method which can 
be obtained from this differential system becomes relatively slow for the zeros of 
Pi"^*^' close to ±1 [4]. Be cause the extreme zeros approach to ±1 as n ^ +oo, the 
efficiency for the computation of such zeros decreases as the order increases. Similar 
problems arise, for example, when a — > — 1+ or /? — > — 1+. In fact, the number of 
iterations required to compute the extreme zeros tend to infinity in these limits. 
Similar problems take place for Laguerre polynomials L'^{x) for the smallest (pos- 
itive) zero. Fortunately we will later show how the selection of Wu{x) as contrast 
function gives a much better asymptotic behavior for the resulting fixed point it- 
eration for the extreme Jacobi zeros. For the Laguerre case, a similar solution is 
possible. 

These two examples illustrate the need to analyze the convergence of the result- 
ing fixed point iteration for the different available contrast functions. Although 
for any adequate contrast function (satisfying the necessary conditions [7, 1]) the 
resulting fixed point method is quadratically convergent, the non-local behavior 
of the method and the corresponding estimation of first guess values for the zeros 
may result in disaster for certain contrast functions in some limits. As a result of 
this study, we will obtain explicit methods for the computation of the real zeros 
of hypergeometric functions with a good asymptotic behavior and a fairly uniform 
convergence rate in the whole range of parameters. 



2. THEORETICAL BACKGROUND 
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Let us now briefly outline the main ingredients of the numerical method. For more 
details we refer to [1, 7]. It was shown in [1, 7] that, given a family of functions 
{yk^\y'k^}, depending on one parameter k, which are independent solutions of 
second order ODEs 

y'k + Bk{x)y'k + Ak{x)yk =0,k = n,n-l (3) 

and satisfy relations of the type: 

yn= an{x)yn + dn{x)yn-i 
y'n-i= bnix)yn-i + en{x)y„ 

the coefficients a„(a;), bn{x), d„{x), e„(a;), Bk{x) and Ak{x) being continuous and 
dn^n < in a given interval [a;i,a;2], fixed point methods (Eq. (7)) can be built to 
compute all the zeros of the solutions of (3) inside this interval. These difference- 
differential equations (4) arc called general because they arc satisfied by a basis 
of solutions {y'k^y^^}- The fact that the DDEs are general and with continuous 
coefficients in an interval / implies [7] that rf„e„ ^ in this interval. Conversely, 
given {yn\yn-i}, i = 1,2 independent solutions of the system (4) and (i„e„ ^ 0, 
then {y'k\y^^}, k = n,n — 1, are independent solutions of the ODEs (3). The 
method can then be applied to compute the zeros of any solution of such ODEs. 
It was shown that the ratios Hi{z) {i = ±1): 

i^.(.)^-.sign(.„jK„.£^ 

/ ^ \V2 (5) 

Kn, = {--^\ , z{x) = J y/-dn,en,dx 

where n+i = n + 1 , n_i = n, satisfy the first order equations 

Hi{z) = 1 + Hi{zf - 2ni{x{z))Hi{z) (6) 

where 

mix) = i I ( a„. -hn, + \ (-^ - ^ 

and the dot means derivative with respect to z while the prime is the derivative 
with respect to x. Using Eq. (6), one can show that 

Ti{z) = z- BxciB.n{Hi{z)) (7) 

are globally convergent fixed point iterations (FPI): given a value Zq between two 
consecutive zeros of yni{x{z)) (consecutive singularities of Hi(z)), the 

iteration of (7) converges to z„ G (z^.,^^.), where a;„ = x{zn) is a zero of yn{x)- 

Global bounds for the distance were provided which lead to iteration steps that 
can be used to compute new starting values for obtaining all the zeros inside a given 
interval. 

It was shown that, in intervals where r/i does not change sign, either 

- Zn\ < 7r/2 and \z{xl^ - z^\ > 7r/2 > 0) (8) 
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or 

\z{xl.) - Zn\ > 'Jt/2 and \z{xl.) - z„\ < n/2 {r]i < 0). (9) 

In this way 7r/2 is the choice for the iteration step when the second situation takes 
place {rii < 0); this means that if yni{x{z)) has at least a zero larger than Zn, then 

lim r(^)(z„ + 7r/2) (10) 

is the smallest zero larger than Zn- Similarly, when r]i > 0, the iteration step will 
be —tt/2 instead of tt/2 (backward sweep). When rji changes sign, forward and 
backward schemes can be combined [1]. 

The functions Hi{z) can be written as a ratio of functions Hi{z) = yn{z)/yni{z), 
where yk{z) = Xk{z)yk{x{z)) and have no zeros, in such a way that the functions 
yn (z) and (z) (with the same zeros as yn{x{z)) and y^ {x{z)) respectively) satisfy 
second order ODEs in normal form; 

(^^Vn _|_ A if — n ^ -L A f/ — n 

kn{z) = l+rii-r]f, kn,{z) = 1 - - r?? . 

Finally, wo recall that using monotony conditions of An{z) the iteration steps 
±7r/2 (Eq. (10)) can be improved according to Theorem 2.4 of [1]: 



Theorem 2.1. // z-\ < zq < zi are three consecutive zeros of yn{x{z)) 
r]i{z)An{z) > in (z_i, zi) then zj = lim„^oo T^"'\zo + Azo) where Azq = zq — Z-j, 
j = sign{r]). The convergence is monotonic. 



2.1. Oscillatory conditions 

Wc arc interested in computing zeros of oscillatory solutions of second order 
ODEs and, in particular, on building algorithms for the computation of the zeros 
of the hypergeometric functions. If a second order differential equation has a given 
number of singular regular points, we divide the real axis in subintervals determined 
by the singularities and search for the zeros in each of these subintervals. We only 
apply the algorithms if it is not disregarded that the function can have two zeros 
at least in the subinterval under consideration. 

We consider that an ODE has oscillatory solutions in one of these subintervals if it 
has solutions with at least two zeros in this subinterval; otherwise, if all the solutions 
have one zero at most we will call these zeros isolated zeros. The fixed point 
methods (FPMs) before described deal with the zeros of any function satisfying 
a given differential equation, no matter what the initial conditions are on this 
function. Isolated zeros for a given solution depend on initial conditions or boundary 
conditions for this solution and are, in any case, easy to locate and compute. 

There are several ways to ensure that a solution yn{x) has at most one zero in 
an interval; among them: 



COMPUTATION OF ZEROS OF HYPERGEOMETRIC FUNCTIONS 



5 



Theorem 2.2. // one of the following conditions is satisfied in an interval I 
(where all the coefficients of the DDEs are continuous) then y„ and (Um) have at 
most one zero in the interval I (trivial solutions excluded): 

1. dn,{x)eni{x) >0 inl [7]. 

2. \r],{x)\ > 1 inl [1]. 
3.1„ < (Ar,, < 0) [IJ. 



The condition Cn-dn^ < is required for the method to apply. Furthermore, it 
is known that when the DDEs (4) are general, d„.e„j can not change sign. There- 
fore dmCrii < is a clear signature for the oscillatory character of the differential 
equation. 

2.2. Hypergeometric functions; selection of the optimal DDEs 

For hypergeometric functions several DDEs are available for the construction of 
fixed point iterations (FPLs), depending on the selection of contrast function. 

Let us start by considering, for example, the case of the confluent hypergeometric 
equation 

xy" + (c - x)y' -ay = (12) 
One of the solutions of this differential equation are Kummer's series 

^k 



M{a, c, x) =iFi (a; c; x) = ^ 



for which different difference-differential relations are available. Indeed, denoting 
an — a + kn, = J + rnn and yn = M(a„,7„,a;) we will have different sets of 
DDES (Eq. (4)) for different selections of {k,m). 

For Gauss hypergeometric functions 2Fi(a, 6; c; a;), which are solutions of the 
ODE 

x{l - x)y" + [c-{a + b+ l)x\y' -aby = 0, (13) 

the possible DDEs are determined by three-vectors with integer components, that 
is, we will consider ?/„ = 2F1 {a + kn, (3 + In;^ + run; x) and the associated DDEs 
will be named {k, I, m)-DDEs. Finally, for the case of the hypergeometric functions 
oFi(;c;x), we can only consider families y„ = oFi(;7 + kn;x) and the different 
relations are described by the integer numbers k. 

Our FPMs can only be applied to solutions of second order ODEs. This restricts 
our study to the hypergeometric functions oFi(; c; x), 2Fo(a, b; ; x), iFi(a; c; x) and 
2Fi{a,b;c;x). 

Regarding the selection of the different DDEs available, we will restrict ourselves 
to: 



1. DDEs with continuous coefficients except at the singular points of the defining 
differential equations. 
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2. The most simple DDEs in a given recurrence direction which allow the use of 
improved iteration steps. Taking as example the case of confluent hypergeometric 
functions, this means that the (1, 0)-DDE will be described, and the analysis of the 
(-1,0), (2,0),... DDEs will be skipped. 

The first restriction is convenient for simplicity and it means that the problem 
function and the contrast function have zeros interlaced in each subdivision of 
the real interval defined by the singular points of the differential equation; this 
is a convenient property for a simple application of the FPMs and enables the 
application of each DDE to compute all the zeros in the different subintervals of 
continuity of the solutions of the differential equation. 

Regarding the second restriction, and considering the case of confiuent hypergeo- 
metric functions as example, it should be noted that for the (fc, m)-DDE, generally 
two FPIs are available, one of them based on the ratio H-i = j/„(x)/y„_i(a;) and 
a second one based on H+i{x) = yn{x)/yn+i{x). As described in [1], generally 
one of these two iterations is preferable because improved iteration steps can be 
considered according to Theorem 2.1 (Theorem 2.4 in [1]). If we considered the 
(— /c, — to)-DDE, the two associated ratios Hi, i = ±1, would be the same as be- 
fore (replacing i by —i). Because both selections of DDEs are equivalent, only one 
of them will be discussed. By convention, we will consider pairs (k, m) for which 
the iteration on can take advantage of the monotony property of A„(.x(z)), 
as described in [1], for classical orthogonal polynomial cases (Jacobi, Hermite, La- 
guerre). Once we have fixed this criterion the index i in Equations (6)-(ll) can 
be dropped. We consider the following additional notation: given a vector il with 
integer components, we will denote by DDE(u) (FP(u)) the corresponding DDE 
(FPI) based on the ratio yn/Vn-i for a; > 0. 

On the other hand, and considering the confluent case as illustration, we will 
not analyze DDE(2,2) nor any successive multiples of the DDE(1,1). This is so 
because the first restriction is generally violated if successive multiples of a DDE 
are considered (there are exceptions to this; see the case of qFi hypergeometric 
functions). 

2.2.1. Selection of the optimal DDEs 

There are several (and related) criteria to select among the available DDEs to 
compute the zeros of a given function y. The associated FPMs tend to be more 
efficient as wc are closer to any of the following two situations: 

1. rj{x) = 0, 

2. the coefficient A„ is constant. 

Of course, the first condition implies the second one (see Eq. (11)). The first 
condition makes the FPI converge with one iteration for any starting value. The 
second condition makes the method an exact one using improved iteration steps 
(there is even no need to iterate the FPM). 

Let us recall that the FPIs associated to a given system of DDEs are quadrati- 
cally convergent to a zero zq (in the transformed variable z) with asymptotic error 
constant ri{x{zo)). Therefore, the smaller |??(x)| is the fastest the convergence is 
expected to be, at least for starting values close enough to zq. On the other hand, 



COMPUTATION OF ZEROS OF HYPERGEOMETRIC FUNCTIONS 



7 



the smaller the absolute value of variation of A„ is, the better the improved iter- 
ation (Theorem 2.1) will work because this implies the exactness of the iteration 
criteria to estimate starting values from previously computed zeros. This second 



however, more relevant to improve the iterative steps for obtaining starting values 
to compute zeros than to improve the local convergence properties of the fixed point 
methods, which are quadratically convergent anyway. 

Indeed, as was described in [4], the natural FPIs for orthogonal polynomials of 
confluent hypergeometric type (FP(— 1, 0)) tend to converge slowly for the compu- 
tation of the first positive zeros when they become very small. This, for instance, 
is the case for Lagucrrc polynomials L"{x) when a — 1+. The reason for this 
behavior lies in the fact that the associated change of variables is singular at a; = 0: 



In this way, the interval of orthogonality for the Laguerre polynomials (0, +oo) 
is transformed into (—00,00) in the z variable. This means that the zeros which 
are very small in the x variable, tend to go to —00 in the z variable. Therefore 
after computing the second smallest zero, X2, the next initial guess for the FPI, 
z{x2) — 7r/2, may lie well far apart for the value z{xi), Xi being the smallest zero. 
Although it is guaranteed that the FPI will converge to z{xi), it could take a 
considerable number of iterations to approach this value. 

Fortunately, we will see that the rest of FPMs (different to FP(fc,0)) do not 
show such a singularity; therefore, we expect better behavior near a; = for these 
iterations. 

This suggests that, given two FPIs with associated change of variables zi{x), 
Z2 {x) respectively, one should choose that one which gives the largest displacement 
in the x variable for the same step in the corresponding z (the typical value being 
7r/2). Let us stress that the possibility of passing the next zero is ruled out: in the 
algorithms the sequence of all z values calculated in a backward (forward) sweep 
form monotonically decreasing (increasing) sequences. 

We will therefore say that the change of variables zi behaves better than Z2 if 
x{zi + Az) > x{z2 + Az), for a typical value of Az (w 7r/2). Given the definition 
of the changes of variables z{x): 



where <i„ = \/—dnen, we can say that the change of variable zi{x) (and its asso- 
ciated FPI) is more appropriate than the change Z2{x) when its coefficient d„ = 
V—dnCn is smaller than the corresponding coefficient for Z2{x). 

Therefore, an alternative non-local prescription to that one dealing with A„ is the 
following: among the possible DDEs and associated fixed point iterations, choose 
that one for which |(i„e„| is smallest. As we will see, this is an easy to apply 
criterion which correctly predicts the more appropriate DDEs and FPI depending 
on the range of the parameters and the dependent variable. 



criterion (on the variation of A) is more difficult to apply, as we will later sec. It is. 



z = 



\/{b— a)(l — a) \nx 




dz 
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3. ANALYSIS OF HYPERGEOMETRIC FUNCTIONS 

We will use the DDEs satisfied by hypergeometric series as generated by the 
Maple package hsum.mpl [5]. The results for each of the family of functions (change 
of variable z(x), function rj(x), etc.) considered can be automatically generated 
using the package zeros. mpl [4]. 

3.1. Hypergeometric function oFi(;c;a;) 

The ODE satisfied by the function y{x) = oFi(; c; x) is 

x'^y" + cxy' -xy = (14) 

The solutions of these differential equations have an infinite number of zeros for 
negative x and are related to Bessel functions: 

oFi(; c; z) = r(c)(-0)(i-^)/Ve_i(2V^) 



3.1.1. First DDE 

Let us consider the DDEs for the family of functions t/„ = oFi(; 7 + n; —x) 
The DDEs for this family reads: 

y'n ^ ^^^-^Un-l + ^-X-^Un ^^^^ 

y'n-1 — ~ Q — \ yn 

where c = 7 + n. 

The relation with Bessel functions can be expressed saying that, if y{x) is a 
solution of x^y" + {u + l)xy' + xy = then 

w{x) = x^yix^/A) (16) 

is a solution of the Bessel equation x^w" + xw' + {x^ — v^)w = 0. 

Transforming the DDEs (15) as described in [1], and in Section 2, the relevant 
functions are: 

r?(^(x)) = A„ = l-^^^^, z{x) = 2^ (17) 

The fixed point method deriving from this set of DDEs will have identical per- 
formance as the system considered in [1], which holds for Ricatti-Bessel functions 
jv{x) = \/xJv{x), that are solutions of the second order ODEs y" + A{x)y = 0, with 
A{x) = 1 — (i^^ — 1 /4) / . The identification of both methods with the replacement 
X x^ /A (according to Eq. (16) and to the change of variable z{x)) is evident 
by comparing the A{x) and A„(a;) coefficients. This is not surprising given that 
both methods compare the same problem function, Ji/{x), with the same contrast 
function, .J^^i(x), up to factors which do not vanish (for example, the factor -^/x 
for Ricatti-Bessel functions) and up to changes of variable. 



3.1.2. Second DDE 
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For this type of hypergeometric functions, the only alternative DDEs that can 
be built are those based on the family of functions t/„ = oFi(; 7 + mn; — x); m = 1 
corresponds to the DDEs (15). We will only consider the case m = 2 (equivalent, 
in the sense described before, to m = —2). For |m| > 2, the DDEs violate the 
first imposed condition on the continuity of the coefficients. For the functions 
Vn = oFi(; 7 + 2n; —x), the associated DDEs read: 

/ _ (c-2)^ + (c-2)-a: , c-1 

{c-2)x ^""^ X Vn-i 

(18) 

y'n-i = --^Vn-i - (c_l)(c_2)^^" 
where c = 7 + 2n. The relevant functions in this case are (again, writing u = c—1): 

il{z{x))=^^^^-l 

K{z{x))={l^)\Ax-{y''-l)) (19) 
z{x) = ^ 



3.1.3. Comparison between DDEs 

The second DDE is no longer equivalent to the first one; in fact, it has quite 
different characteristics. Let us compare the expected performance of these two 
DDEs, according to the different criteria described above. 

To begin with, the r]{x) parameter never vanishes for the first DDE (DDEl), 
except when v = 1/2, in which case the method with improved steps is exact 
without the need to iterate the FPI even once (forward or backward sweeps [7] are 
used depending on the sign of r]{x)). In contrast, DDE2 has an 77-function which 
changes sign at Xn = {v — 1)^/2 and the zeros are computed by an expansive sweep 
[7]. Close to Xff wc can expect that DDE2 tends to behave better in relation to 
local convergence, because the asymptotic error constant tends to be small. 

We observe that, as x ^ +00, r]{x) goes to zero for DDEl but it tends to — 1 
for DDE2. This suggests that DDEl will have faster local convergence than DDE2 
for large x. On the other hand, as v increases, ri{x) becomes larger, however, it is 
difficult to quantify the impact on local converge because as v increases also the 
smallest zero becomes larger. Let us also take into account that DDE2 will have 
small ri{x) for x close to x^, which becomes large for large v. 

Regarding the behavior of An{z), it is monotonic for DDEl and has a maximum 
for DDE2, which allows the use of improved iteration steps. For DDEl, An{z) is 
constant when u = 1/2, which means that sweep with improved iteration steps is 
exact, as commented before. The maximum for DDE2 is at Xm = {v^ — l)/2, where 
An{xm) = ~ l)/(i^ + !)• Around this extremum, the improved iteration steps 
tend to work better because A„(a;) will be approximately constant; how constant 
A„(x) is around Xm can be measured by the convexity at this point. We find: 
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where Zm — z{xm)- As ly becomes larger, A„(z) becomes smaller around the maxi- 
mum of A„(z) and the improved iteration will work better. This fact again, favors 
DDE2 for large v. 

Finally, considering the criterion of smaller Z)„ = |rf„e„|, we find that, for DDEl 

Dn = - (20) 

X 

while for DDE2 

Dn = ^ (21) 

This again shows that DDEl will improve as x increases while DDE2 will be 
better for largo z/. Numerical experiments show that for v > 100 the second DDE 
is preferable over the first, particularly for computing the smallest zeros. 

The different criteria yield basically the same information. However the prescrip- 
tion on Dn = |rfne„| is the simplest one to apply. Prom now on, we will not repeat 
the analysis for the different criteria. Instead, we adopt this last criterion to analyze 
the rest of cases. 




600 800 1000 



Figure 1. Left: Ratio between the number of iterations needed for the second and 
first DDEs for the computation of the zeros of oFi(; 11, —x) (the zeros of Jio{2^/x)). 
Right: Number of iterations needed for the first DDE. 
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Figure 2. Left: Ratio between the number of iterations needed for the second 
and first DDEs for the computation of the zeros of oFi(; 201, — a;) (the zeros of 
•/20o(2\/^))- Right: Number of iterations needed for the first DDE. 

3.2. Confluent hypergeometric function 

For the confluent hypergeometric case we have a larger variety of DDEs to choose, 
because we can choose families of functions j/„ = iFi{a+ kn;c + mn;x), or, more 
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generally, y„ = (j){a + kn;c + mn;x), being <p any solution of Eq. (12). The families 
which give rise to DDEs satisfying all our requirements are three, corresponding to 
the following selections of {k, m): (1, 0), (1, 1) and (0,-1). 

We will give the corresponding DDEs and the associated functions. At the same 
time, we will restrict the range of parameters for which the functions are oscillatory 
(considering the first oscillatory condition in Theorem 2.2). 

We can restrict the study to a; > because, if (/)(a; c; x) is a solution of Eq. (12), 
then e'^(/>(c — a; c; —x) is also a solution of Eq. (12). 

3.2.1. {k, m) = (1, 0) -> y„ = (pfa + n; 7 ; x) 

Let us write, for shortness and in order to compare with other recurrences a = 
a + n, c = -y. As before commented, we consider simultaneously the equivalent 

directions (k,ni) = (1,0) and {k,m) — (—1,0) but we present only the DDEs and 
related functions for the recurrence direction for which the FPI based on the ratio 
H-i = Dn/yn-i can be used with improved iteration steps (Theorem 2.1); this is 
the direction (1,0). 

The (—1,0) direction is the natural one for orthogonal polynomials of hypcrge- 
ometric type (Laguerre, Hcrmite), which are related to confluent hypergeometric 
series of the type iFi(— n;7;a;) (see, for instance, [6], Eqs. (9.13.8-10)). 

The DDEs for (fc,m) = (1,0) read 



/ ^ a-c+x _ a-c 

Un X X yn-l 



yn-1 — X i*"-! ^ X 



(22) 



and therefore, applying the first oscillatory condition of Theorem 2.2, the parame- 
ters are restricted to: 

(a - c)(a - 1) > 

otherwise the functions will be non-oscillatory. If we repeat the same 

argument for (k,m,) = (—1,0) or, equivalently, apply the same criteria for the 
DDEs for the functions y„, yn+i the following restriction is obtained: 

(a - c + l)a > 0. 

The associated functions for this DDE are 

^(^(a;)) 2a + c + x 



2^/{c-a){l-a) 

An{z{x)) = 



x^ + 2{c - 2a}.v - [c- l f (23) 



4(c-o)(l-a) 



z{x) = i/(c — o)(l — a) Inx 

Let us notice that ■r]{x) becomes negative for large x, which gives the most ap- 
propriate sweep for large x (forward) since A„(z(x)) decreases for large x. 

Observe that the lack of a singularity in An{z{x)) is only apparent because the 
function z{x) is singular at a; = 0. 
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3.2.2. {k, m) = (0, -1) y„ = (p(a; j - n ; x) 

As before, we denote a = a, c = ^ — n. The DDEs for the system read: 

y'n= yn + ^^Vn-l 

(24) 

= — + xyn 

which imphes the restriction {c—a)x > for the solutions to have oscillatory nature. 
Repeating the same for {k, m) = (0, 1), we arrive at the condition (c — a — \)x > 0, 
which for x > gives a more restrictive condition c — a > 1 (for x < 0, the first 
condition is more restrictive and gives a — c < 0). 

Let us recall that the condition c— a>lforx>0 means that, if this condition 
is not met, neither j/„ nor y„_i can have two zeros in a; > 0. Given that we are 
interested in computing zeros of oscillatory functions, we consider these restrictions. 
Let us however notice that the possible isolated zero of y„ for x > and < c — a < 
1, could be also computed by means of the FPM associated to the DDEs (24). 

Considering also the restrictions imposed in the previous selection of DDEs, we 
obtain the following: 

Theorem 3.1. Let y be a solution of the confluent hypergeometric equation (12) 
for x> Q. If y has at least two zeros then c — a > 1 and a < 0. 

Let y he a solution of the confluent hypergeometric equation (12) for a; < 0. If y 
has at least two zeros then c — a < and a > 1 . 

For more detailed results on the number of zeros of confluent hypergeometric 
functions, we refer the reader to [2], Volume 1, Section 6.16. 
The associated functions are 

4-\/(c — a)x 

~K{z{x)) = 8c.r-lGxa ,3 + 8c-4c^-4x^ (25) 
^ ^ " 16(c — a)x 

z = 2i/ (c — a)x 

Let us notice that r]{z{x)) becomes negative for large x, which gives the most 
appropriate sweep for large x (forward) since, for positive x, An{z{x)) decreases for 
large x. 

3.2.3. {k,m) = (1, 1) ^ j/„ = (l)(a + n; 7 + n ; x) 

The DDEs are the following (writing a = a + n, c = 7 + n). 

2/n = ^ X ~y^ ~^ yn-1 

(26) 

yn-1 — Q — \yn 

which implies the oscillatory condition (a — l)x > 0; considering {k, m) = (—1, —1) 
or, equivalently, the DDEs relating y„, yn+\ and their derivatives, the condition 
obtained is aa; > 0. These conditions are consistent with Theorem 3.1. 
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The functions associated to these DDEs for a; > are 



V 



2x + 3- 2c 
V(l-a)a; 



z 



2^{l-a)x 



(27) 



An{z{x)) = 



16xa + 4a:^ - 8xc + 3 - 8c + 4c^ 
16(-l + a)a; 



This iteration cannot be used for c = 1 (sec Eq. (26)) unless the FPI stemming 
from the ratio yn/Un+i is used, in which case the improved iteration steps cannot 
be apphed. 

3.2.4- Comparing fixed point iterations 

As commented before, we will consider the prescription consisting in choosing 
the DDEs for which the product D„ = — (i„e„ is smaller. Applying literally this 
criterion, we find the following preferred regions of application; 

1. FP(1,1) should be applied for a; < c — a and FP(1,0) for a; > c — a. 

2. FP(1,1) is always better than FP(0,-1). 

3. FP(0,-1) is better than FP(1,0) for x < 1 - a, but 1 - a < c - a and FP(1,1) 
is better there. 

Therefore, the best combination of the considered FPIs is FP(1,1) for x < c — a 
and FP(1,0) for a; > c — a. The iteration FP(0,-1) has the same behavior as 
FP(1,1) and can be used as a replacement when c = 1 (in this case FP(1,1) can 
not be used). In fact, the FP(0,-1) and FP(1,1) are not independent because, as it 
is well known, if 'tp{a;'y;x) are solutions of the confluent hypergeometric equation 
xip" + (7 — x)^p' — atl) = 0, then y{a; c; a;) = x^~'^il){l + a + c; 2 — c; x) is a solution 
of xy" + (c — x)y' — ay = 0. 

Let us illustrate this behavior with numerical examples: 

In Fig. 3 we compare FP(1,1) with FP(1,0) for the case of Laguerre polynomials 
L"(a;) for the quite extreme case n = 50, a = —0.9999. On the left, the ratio 
between the number of iterations employed by FP(1,1) and FP(1,0) is shown as 
a function of the location of the zeros. The first two zeros are skipped in the left 
figure to show in more detail the behavior for large x (for the first zero the ratio was 
40 while for the second it was 5). The improvement for small x when considering 
FP(1,1) is quite noticeable. For larger x {x > c — a:^ 50), FP(1,0) works generally 
better than FP(1,1) but the improvement is not so noticeable. On the right figure, 
the number of iterations used to compute each zero when considering FP(1,1) is 
shown. 
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Figure 3. Left: Ratio between the number of iterations needed by FP(1,0) 
and FP(1,1) for the calculation of the zeros of the generalized Laguerre polynomial 
L^Q^^^^{x), as a function of the location of the zeros. Right: Number of iterations 
needed by FP(1,1). 

In conclusion, FP(1,1) is a more appropriate choice than the natural recurrence 
for orthogonal polynomials (FP(1,0)), although FP(1,0) sUghtly improves the con- 
vergence of FP(1,1) for a; > c — a. 

In Fig. 4 we compare FP(1,0) with FP(0,-1) for the case of generalized Laguerre 
polynomials but now for a choice of the parameters n = 50, a = 0. This situation 
corresponds to the case c = 1 where FP(1,1) can not be applied. As in Fig. 3, the 
ratio of the number of iterations for the first two zeros is not plotted (the ratio of 
the first zero was 8 and for the second it was 5). As expected from our previous 
analysis, the iteration FP(0,-1) behaves quite better than FP(1,0) for small x. 
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Figure 4. Left: Ratio between the number of iterations needed by FP(1,0) and 
FP(0,-1) for the calculation of the zeros of the generalized Laguerre polynomial 
L5o(x), as a function of the location of the zeros. Right: Number of iterations 
needed by FP(0,-1). 



3.3. Mysterious hypergeometric function 

This is the name given to the hypergeometric scries 2Fo(a, b;;x), which diverges 
for all a; 7^ (except in the terminating cases) and can only be interpreted in an 
asymptotic sense. It is well know that, for negative x we have: 



2Fo(a, 6; ; a;) = (-l/a;)"f/(a, 1 + a-b, -l/x) 



(28) 
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where U{a,c,x) is a solution of the confluent hypergeometric equation (12). The 
functions 2F0 can be analytically continued to the whole complex plane cut along 
the line 3?(z) > 1, = 0. The function (28) is a solution of the 2-0 hypergeo- 
metric differential equation: 

x^y" + [-1 + x{a + b+ l)]y' + aby = 0. (29) 

In general terms, without referring to any particular solution of the correspond- 
ing differential equations, the problem of computing the real zeros of the mysterious 
hypergeometric function for negative x can be transformed into a problem of com- 
putation of the zeros of confluent hypergeometric functions Eq. (28). This is so 
because one can check that, if we denote by y{a, (3, x) a set of solutions of the con- 
fluent hypergeometric equation, then iu{x) = \x\~'^y{a, 1 + a — b,—l/x), for x > 
or a; < 0, arc solutions of Eq. (29). 

For this reason and for brevity we omit further details. 

3.4. Gauss Hypergeometric Functions 

Let us consider the hypergeometric fimction 2Fi(a, 6; c; x). We will consider 
the DDEs for families of functions of the type ip{a + kn, /? + mn; 7 + In; x), with 
ip{a,b;c;x) solutions of the hypergeometric equation (13); we use the DDEs for 
2F1 (a, &; c; x) series as generated by hsum.mpl. 

Similarly as we did for the confluent case, we can obtain oscillatory conditions 
for the coefficients a, b and c, depending on the range of x. If these conditions are 
not satisfied by the parameters, then we can assure that if there exists one zero of 
the fimction, this is an isolated zero. As before, these conditions are obtained by 
requiring that d„e„ < 0; combining the restrictions imposed by this condition for 
the DDEs that we will later show we obtain the following: 

Theorem 3.2. Let tp{a, b; c; x) be a solution of Eq. (13) defined in (—00, 0), then 
if this function is oscillatory in this interval (it has at least two zeros) then one of 
the following sets of conditions must be verified: 

( (CI): a<0,b<0,c-a>l,c-b>l 
\ (C2): a>l,b>l,c-a<Q,c-b<Q 

Similarly, if tp{a,b;c;x) is a solution of Eq. (13) defined in (0,1), the oscillatory 
conditions are 

( (C3): a<0,6>l,c-a>l,c-6<0 
\ (C4): a>l,b<Q,c-a<Q,c-b>l 

Finally, if 4'{o,,b: c; x) is a, solution of Eq. (13) defined in (l,+oo) and is oscilla- 
tory in this interval, one of the following sets of conditions must be verified: 



( (C5): a<0,6<0,c-a<0,c-6<0 
\ (C6): a>l,6>l,c-a>l,c-6>l 
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The different sets of conditions for the three subintervals in Theorem 3.2 can 
be obtained by combining the restrictions obtained from the following values of 
{k,l,m): (±1,0,0), (±1,±1,0), (±1,0, ±1), (±L=f1,0), (0,0, ±1) and (±1, ±1, ±1). 
In any case, the analysis for the three different subintervals are not independent 
because, as it is well know (see [2], Vol. I, Chap. II), if we denote by tp{a, /?; 7, x) the 
solutions of the hypergeometric equation x{l — x)y" + {j — {a+ l3+l)x)y' — a(3y = 
in the interval (0, 1) one can write solutions in the other two intervals by using that 
both 

y{a, b; c; a;) = (1 - x)~°'ij){a, c-b;c; x/{x - 1)) , a; < (30) 

and 

y{a, b; c; x) = x~'^'ip(a, a+1 — c;a + 6+ l — c;l — 1/x) , a; > 1 (31) 

are solutions of the hypergeometric differential equation x(l — x)y" + (c — (a + & + 
l)x)y' — aby = 0. With this, it is easy to see that the conditions CI and C2 can be 
obtained from C3 and C4 respectively by using Eq. (30), while C5 and C6 derive 
from C3, C4 and Eq. (31). 

Notice, in addition, that the 6 conditions in Theorem 3.2 are mutually exclusive, 
which means that 

Theorem 3.3. Given three values of the parameters a, b and c, at most one of 
the subintervals (— oo,0), (0,1), (l,+oo) possesses oscillatory solutions (solutions 
with at least two zeros). 

The oscillatory conditions C6 is of no use for non-terminating hypergeometric 
series, because they diverge for a; > 1; however, this is a possible case for other 
solutions of the differential equation. The conditions C3 correspond to Jacobi poly- 
nomials Pj,"''^\x) = 2Fi(-n, 1 + « + /? + ri;a + l:(l-. t)/2) (a,/3 > -1) 
of order n > 2 (for order n = 1 the conditions are not satisfied, not surprisingly 
because according to our criteria a polynomial of degree 1 is non-oscillating). Par- 
ticular cases of Jacobi polynomials are Gegenbauer {a = (3), Legendre (a = /? = 0) 
and Chebyshev (a = /3 = — 1/2) polynomials. 

3.4-1- DDEs and change of variables 

In this section, we compile the expressions for the different DDEs as well as the 
associated change of variable. For brevity, the associated functions r]{x) and A„(x) 
are not shown. 

It is understood that a = a + kn, b = f3 + In, c = 'y + mn for DDE(fc, /, m). 
1. DDE(1,0,0) 

1,' — -a-bx + c , a - c 

- x{x-l) + x(x-l)2/"-i 

11' — 1 ~ 0" !i , 4- 0'~ ^ 11 
Vn-l — X yn-l T 2; (/n 

Change of variable: z{x) = —2^{c — o)(l — a) tanh~^(v'l — x) 
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2. DDE{1,1,Q) : 

, _ / -{a + b)(a + b-c-l) + ab-c g + b-A „ , (a-c)(c-b) 

x{a + b-c-l) ^ a;(l-a;)^^"^a;(l-x)(a + 6-c-l)^"-^ 

/ _ a + b-ab-1 „ (l-x)(a- 1)(1 -b) 

^"-1 ~ (a + 6-c-l)a;^"-^ (a + 6 - c - l)a; ^" 



Change of variable: = - c)(c - a)(6 -1)(1 - a) ^ 
^ ^ ' \a + b — c — 1\ 

3. DDE{1,1,2) : 

, _ (l-g;) [(i-a-b){c-l) + ab] + {c-l){l + a + b-c)-ab _ (1 - c) 

a;(l-x)(c-2) ^" a;(l - x) ^/"-i 

,/ _ l-a-6 + a6 „ , _ a;(a - c + 1)(1 - a)(c - b - 1)(6 - 1) 
^"-1- (l-x){c-2)^"-i (l-:r)(c-l)(c-2)^ ^" 



Change of variable: z{x) = _^ {c - a - m -a)il + b - c){b - 1) _ 
4. £)£)i;(l,0,l) : 

/ _ 1 + xb - c _ (J- - c) „ , 

, _ (l-g) (i-a)(6+l-c) 
y^-i- (1-a;)^"-^^ (l-a;)(l-c) ^" 



Change of variable: = 2-^/ (1 — o)(6 + 1 — c) tanh ^ y^. 

5. D£>i;(l,-1,0) : 

/ _ j^ x{b — a + l) + a — c b{c — a) 

^" ~ x{l ~ x){b - a + I) + a;(l-x)(6-a+l)^"-i 

,/ - fl + (l-a)(l + 5-c) 

yn-i- U «J x(l-2;)(6-a+l) ^"-1 a;(l-x)(6-a + l)^" 

Change of variable: z{x) = s/b{c- a){l-^a){l + b - c) 

6. DDE{Q,Q,-l) : 

/ _ 6 + a-c _ (6-c)(c-a) 
2^"- 1-a; (l-a;)c ^""^ 

/ c c 



Change of variable: z{x) = -\/(6^^c)(c^^a) arcsin(2x — 1). 
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7. DDE{1,1,1) : 



Vn = 



_ x{a + b-l)-c+l 



1-c 



(1 - x) ^" a;(l-a;)^"-^ 



, _ (6-l)(l-a) 



Change of variable: z{x) = \/ {b — 1)(1 — a) arcsin(2a; — 1). 



3.4-2. Comparison of FPIs 

The following tabic shows the \Dn\ = \ — dnCn] coefficient for the FPIs (1,0,0), 
(1,1,0), (1,1,2), (1,0,1), (1,-1,0), (0,0,-1) and (1,1,1). 



Iteration 



\D„ 



(1,0,0) 




(a — c)(a — 1) 
x^{l-x) 




(1,1,0) 




{b - c){a - c){b - l){a - I) 




x^{a + b- c-lf 


(1,1,2) 


(c - a - 1)(1 - a)(l + b- c){b - 1) 


{l-xf{c-2f 



(1,0,1) 



(a - l)(l + 6-c) 



x{l - xY 



(1,-1,0) 



b{c-a){a- l)(c-6- 1) 



x\l-xY{a^b~lY 



(0,0,-1) 



(6 — c){a — c) 



x{l — x) 



(1,1,1) 



(6-l)(a-l) 



x{l — x) 



As can be inferred from the table, the most appropriate FPIs in the interval 
(0, 1) are the (0, 0, —1) and (1,1,1) iterations. As commented, some hypergeometric 
functions in this interval with particular values of their parameters are orthogonal 
polynomials (Jacobi and derived polynomials). In the case of orthogonal polynomi- 
als, the "natural" iteration to be considered is (1, —1, 0) (or (—1, 1, 0) equivalently) 
which is not the optimal iteration. In order to illustrate this fact, let us consider 
the evaluation of the zeros of the hypergeometric function 2Fi(— 50, 54; 5/2; a;) in 
the interval (0,1). The zeros of this function correspond to the zeros of the Jacobi 

f3/2 3/2) 

polynomial P^^' '{1 - 2x). In Fi gure 5, we show the ratio between the number 
of iterations needed by FP(1, — 1, 0) and FP(1, 1, 1) as a function of the location of 
the zeros. 

On the contrary, the most appropriate iteration in the (l,oo) interval is the 
(1,0,0) iteration. This could also be understood taking into account that the best 
iterations in the interval (0, 1) are FPI(0,0,-1) and FPI(1,1,1) and using Eq. (31). In 
Figure 6 we show the ratio between the number of iterations needed by FP(1, 1,1) 
andFP(l,0,0). 
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The main conclusion for 2F1 hypergeometric functions is that the iteration (1, 1, 1) 
is the preferred one and that (0,0,-1) can be considered as a replacement, with 
similar performance. For the other two intervals, the relations (30) and (31) indi- 
cate that the most appropriate iterations will be (1,0,1) for (— oo,0) and (1,0,0) 
for (1, +00); in any case, the solutions in these intervals can be related to solutions 
in (0,1). 
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Figure 5. Ratio between the number of iterations needed by FP(1,-1,0) and 
FP(1,1,1) for the calculation of the zeros of the hypergeometric function 2Fi(— 50, 54; 5/2; x), 
as a function of the location of the zeros. 
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Figure 6. Ratio between the number of iterations needed by FP(1,1,1) and 
FP(1,0,0) for the calculation of the zeros of the hypergeometric function 2Fi(— 30, —32; —70; x), 
as a function of the location of the zeros. 



4. CONCLUSIONS 

We have developed a detailed study of the performance of the available fixed point 
methods for hypergeometric functions. This will allow the construction of efficient 
algorithms for the computation of the real zeros of hypergeometric functions with 
a good asymptotic behavior. The next table summarized the main results. 
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Function 


Contrast function 


Range of application 


Change of variables z{x) 


oFi(;c; -x) 


oFi(;c_; -x) 
oFi(;c_;-x) 


r c < 100 

\ c> 100 , x > cV2 
c> 100, a; < 


2^ 
x/{^-l) 


iFi{a;c;x) 


iFi(a_;c_;a;) 
iFi (o_;c; x) 


x < c — a 
x > c — a 


2^{l-a)x 
-/Vac In a; 


2Fi{a,b;c;x) 


2Fi{a-,b-;c-;x) 




Nbc arcsin(2a; — 1) 



where Nac = a/(c — o)(l — a) and N^c = \/{h— 1)(1 — a) and o_ = a — 1, 6_ = 

6 — 1, c_ = c — 1. 

In the table we restrict ourselves to a; > in all cases, with the additional 
restriction x < 1 for the 2F1 functions. For the rest of the intervals, as commented 
before, relations are available which map these other regions into the intervals 
considered in the table. 

It should be noted that, when a priori approximations to the roots are available 
(for instance, asymptotic approximations like in [8]) the performance of the algo- 
rithms can be improved. However, the methods presented here have the advantage 
of being efficient methods that do not require specific approximations for specific 
functions (which, on the other hand, are difHcult to obtain for three parameter 
functions like the 2F1 hypergeometric functions). In addition, even in the simple 
cases of one parameter functions, the methods are very efficient by themselves. 

To conclude, it is worth mentioning that one on the main reasons for the good 
performance of the algorithms if that the analytical transformations of the DDEs, 
and in particular, the associated change of variable z{x) = J y/—dnendx, tend to 
uniformize the distance between zeros. Generally speaking, the most successful 
methods are those which produce smaller variations of the distances between zeros, 
because the first guesses for the zeros become more accurate. In connection to this, 
these changes of variable lead to interesting analytical information about these zeros 
[3]. 

APPENDIX: MAPLE CODE 

In this appendix we would like to explain how we received the DDEs automat- 
ically. Our Maple code uses a Zeilberger type approach [9] and is completely on 
the lines of [5]. The authors provide a Maple program rules. mpl which can be 
used in combination with hsum.mpl [5] to get equations (15), (18), (22), (24), 
(26), (28), and the DDEs in § 3.4.1. These computations are collected in the 
Maple worksheet rules. mpl. All these files can be obtained from the web site 
http : / /www . mathematik . uni-kassel . de/~koepf /Publikationen. 
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